'''
Created on Jan 30, 2010

@author: oabalbin
'''
import sys
import os
import numpy as np
import scipy as sp
from scipy import stats
from collections import deque

#from numpy.random import permutation


import signatures.parsers.read_gene_lists as gpm
import signatures.stats.Rstats as rstat
from signatures.common.classes import arrayMatrix, cpmMatrix 


def filter_out_constant_genes(expMat, lvar, lmed, log=True):
    """
    returns a list of gen probes with std >= lvar and med >=lmed in the array
    expMat is an object of the class arrayMatrix
    Note: recall this analysis is made at the probe level, and some
    genes have more than one probe
    the expMat.expVal array is a masked array
    """
    
    gen_std = []
    gen_med = []
        
    if log:
        gen_std = np.ma.std(expMat.expVal, axis=1)
        gen_med = np.ma.median(expMat.expVal, axis=1)
    else:
        gen_std = np.ma.std(np.log2(expMat.expVal), axis=1)
        gen_med = np.ma.median(np.log2(expMat.expVal), axis=1)
    
    ng = expMat.expVal.shape[0]
    genes = np.arange(ng)
    var_cut = sp.stats.scoreatpercentile(gen_std,lvar) 
    med_cut = sp.stats.scoreatpercentile(gen_med,lmed)
    
    #print 'ERG sd in samples and sd_cutoff ' + str(gen_std[0]) +'\t'+ str(var_cut)
    gtk_std = set(genes[gen_std > var_cut])
    gtk_med = set(genes[gen_med > med_cut])
    
    gentokeep = np.array(list(gtk_std.intersection(gtk_med)))
    #gentokeep = np.array(list(gtk_std))
    
    #core = ['XRCC5','PRKDC','XRCC6','DDX5','PARP1','AR','NUMA1','PRPF8','RPA1','RPA2','TOP2B']
       
    return gentokeep



def filter_out_constant_and_masked_genes(expMat, lvar, lmed, lmasked, log=True, fold_change=0, tumor_indices=[]):
    """
    returns a list of gen probes with std >= lvar and med >=lmed in the array
    expMat is an object of the class arrayMatrix
    Note: recall this analysis is made at the probe level, and some
    genes have more than one probe
    the expMat.expVal array is a masked array
    """
    
    gen_std = []
    gen_med = []
        
    if log:
        gen_std = np.ma.std(expMat.expVal, axis=1)
        gen_med = np.ma.median(expMat.expVal, axis=1)
        gen_mask = np.ma.count_masked(expMat.expVal, axis=1)
        gen_mean = np.ma.mean(expMat.expVal, axis=1) 
    else:
        gen_std = np.ma.std(np.log2(expMat.expVal), axis=1)
        gen_med = np.ma.median(np.log2(expMat.expVal), axis=1)
        gen_mask = np.ma.count_masked(expMat.expVal, axis=1)

    
    ng = expMat.expVal.shape[0]
    genes = np.arange(ng)
    var_cut = sp.stats.scoreatpercentile(gen_std,lvar) 
    med_cut = sp.stats.scoreatpercentile(gen_med,lmed)
    
    #print 'ERG sd in samples and sd_cutoff ' + str(gen_std[0]) +'\t'+ str(var_cut)
    gtk_std = set(genes[gen_std > var_cut])
    gtk_med = set(genes[abs(gen_med) > med_cut])
    gtk_nomaksed = set(genes[gen_mask < lmasked])  
    
    print 'Number of masked elements '
    print gen_mask
    print lmasked
    print len(gtk_nomaksed)

    print 'Number of median below elements '
    print med_cut
    print len(gtk_med)

    gentokeep0 = gtk_std.intersection(gtk_med)
    
    if fold_change != 0:
        # Change this because I want to specify that it is the mean of a particular group, for example tumors
        # given than the data was already normalized against normals.
        gen_mean_tumor = np.ma.mean(expMat.expVal[:,tumor_indices], axis=1)
        gtk_mean = set(genes[abs(gen_mean_tumor) > fold_change])
        gentokeep0 = gtk_mean.intersection(gentokeep0)
    
    gentokeep = np.array(list(gentokeep0.intersection(gtk_nomaksed)))
    #gentokeep = np.array(list(gtk_std))
           
    return gentokeep
    

def filterout_genesLowLogFoldChange(expMat, fold_change=0, tumor_indices=[]):
    
    ng = expMat.expVal.shape[0]
    genes = np.arange(ng)

    if fold_change != 0:
        # Change this because I want to specify that it is the mean of a particular group, for example tumors
        # given than the data was already normalized against normals.
        gen_mean_tumor = np.ma.mean(expMat.expVal[:,tumor_indices], axis=1)
        gtk_mean = set(genes[abs(gen_mean_tumor) > fold_change])
        gentokeep = gtk_mean.intersection(set(genes))
               
    return gentokeep



def get_variable_genes(expMat, vargenes):
    """
    It needs a expMat object and an list of indicators of gene of interest.
    It returns a list of tuples with index and gene and a gen dictionary
    """
    gn_ind=deque([])
    gen_dic={}
    #expMat.geneSym.pop(0)
    genelist = np.array(list(expMat.geneSym))
    genelist = genelist[vargenes]
    #vg_expMat = expMat.expVal[vargenes,:]
    
    for g, i in zip(genelist,vargenes):
        gn_ind.append((g, i))
        gen_dic[g] = i
    #print gn_ind
    return gn_ind, gen_dic


def get_filtered_expMat(expMat,geneDict,nsamples):
    """
    This functions returns an expression matrix array that was filtered
    according to the genes present in geneDict. This genes usually are 
    the ones that are going to be used in the bfrm analysis.
    """
    
    ngenes = len(geneDict)
    expmat_bfrm = arrayMatrix(ngenes, nsamples)
    
    expmat_bfrm.sample = expMat.sample
    #i=0
    for i, p in enumerate(geneDict):
    #for g, ind in geneDict.iteritems():
        ind = p[1]
        g = p[0]
        
        expmat_bfrm.probeID = expMat.probeID[ind]
        expmat_bfrm.geneSym.append(g)
        expmat_bfrm.geneSynd.append((i,g))  #add((i,g))
        expmat_bfrm.expVal[i,:] = expMat.expVal[ind,:]
        expmat_bfrm.stdVal[i,:] = expMat.stdVal[ind,:]
        #i+=1
        
    return expmat_bfrm

def get_reordered_geneList(mastergeelist,seedgenelist):
    """
    It reorders the rows of expMat, so the n seed genes in seedDict correspond to the first n rows of the matrix.
    This is done for a later use of BFRM 
    Sorting the master list could be eliminated It does not affect anything downstream
    """
    odgenelist = deque([])
    mastergeelist.sort()
    
    ## Check that all elements of seedgenelist are in masterlist
    seedgenelist, mastergeelist = set(seedgenelist), set(mastergeelist)
    seedgenelist.intersection_update(mastergeelist)
    seedgenelist, mastergeelist = list(seedgenelist), list(mastergeelist)
    
    for g in seedgenelist:
        odgenelist.append(g)
        mastergeelist.remove(g)
    
    [odgenelist.append(g) for g in mastergeelist]
        
    return odgenelist
    
    
def get_reordered_samplelist(samplelist,sampleclasses):
    """
    sample_class contains a subset of the samplelist, and each part is classified according to a particular criteria.
    """
    odsamplelist = deque([])
    allsp = set(samplelist)
    
    for spclass in sampleclasses:
        [odsamplelist.append(sp) for sp in spclass]
    
    allsp.difference_update(set(odsamplelist))
    allsp = list(allsp)
    allsp.sort()
    [odsamplelist.append(sp) for sp in allsp]
    
    return odsamplelist
        


def get_matrix_data(inputfile):
    
    matdata = np.loadtxt(inputfile,dtype=float, delimiter='\t', skiprows=1)
    return matdata


def get_header_mat(inputfile):
    
    line = inputfile.readline()
    line = line.strip('\n')
    fields = line.split('\t')
    #print fields
    return fields


def pickledump_matrix_data(expMat, matfile):
    """
    dump a masked matrix array
    """
    np.ma.dump(expMat, matfile)


def pickleload_matrix_data(inputfile):
    """
    load a masked matrix array from dump file
    """
    matdata = np.ma.load(inputfile)
    return matdata


def calc_correlation_mats_deprecated(mat1,mat2, cmpname, corrtype, pvalth, corth):
    """
    mat1, mat, are numpy masked arrays, with the same number of columns
    mat1=compounds
    mat2=pathways
    but the different (equal) number of rows, mat1 > mat2.
    Columns correspond to the same variable, ex: sample
    kend, pear: kendaltau or pearson correlation
    """
    m1 = mat1.shape[0]
    m2 = mat2.shape[0]
    parametric = {'Kend':False,'Pear':True}[corrtype]
    outa_cor=np.zeros([1,6])
    
    npairs=0
    
    nper=int(1000)
    nlabels = int(mat2.shape[1])
    permtlist = array_permutation(nper,nlabels,arr1=[])
    
    for i in range(m1):
        print i
        for k in range(m2):
            outa = np.zeros([1,6])
            GI50 = mat1[i,:]
            path = mat2[k,:]
            if i < m1:
                pcor, pval = calc_corelation_arr(GI50, path, parametric)
                #print 'These are the correlation and pval '+ str(pcor) + '  ' + str(pval)
                
            else:
                break
            
            if pval <= pvalth and corth <= abs(pcor):
                print i, cmpname[i], k, pcor, pval
                
                bpval = corr_permutation_test(GI50, path, permtlist, parametric, pcor)
                
                
                if bpval <= pvalth:
                    print i, cmpname[i], k, pcor, pval, bpval
                    outa[0,0],outa[0,1] = i, k
                    outa[0,2],outa[0,3],outa[0,4], outa[0,5] = pval, pcor, bpval, int(cmpname[i])

                outa_cor = np.vstack((outa_cor,outa))
                npairs+=1
    
    print 'The number of accepted correlation is '+ str(npairs)
    
    return outa_cor[1:,]
    

def calc_correlation_mats(mat1,mat2, cmpname, corrtype, pvalth, corth, permutationsFile):
    """
    mat1, mat, are numpy masked arrays, with the same number of columns
    mat1=compounds
    mat2=pathways
    but the different (equal) number of rows, mat1 > mat2.
    Columns correspond to the same variable, ex: sample
    kend, pear: kendaltau or pearson correlation
    """
    m1 = mat1.shape[0]
    m2 = mat2.shape[0]
    parametric = {'Kend':False,'Pear':True}[corrtype]
    outa_cor=np.empty([1,6])
    
    
    
    nper=int(1e6)
    nlabels = int(mat2.shape[1])
    # permuted labels for the GI50 matrix
    #nlabels = int(mat1.shape[1])
    #permtlist = array_permutation(nper,nlabels,arr1=[])
    if not os.path.isfile(permutationsFile):
        print 'Calculating the permutation matrix'
        permut_matrix= corr_permutation_distr(nlabels, mat1, mat2, nper, parametric=True)
        pickledump_matrix_data(permut_matrix, permutationsFile)
    else:
        print 'Loading the permutation matrix ' + permutationsFile
        permut_matrix=pickleload_matrix_data(permutationsFile)
        
    totc=m1*m2
    this_iter=0
    npairs=0
    for i in range(m1):
        nbe =totc-this_iter
        print nbe
        for k in range(m2):
            outa = np.empty([1,6])
            GI50 = mat1[i,:]
            path = mat2[k,:]
                        
            if i < m1:
                pcor, pval = calc_corelation_arr(GI50, path, parametric)
                #print 'These are the correlation and pval '+ str(pcor) + '  ' + str(pval)
                
            else:
                break
            
            if pval <= pvalth and corth <= abs(pcor):
                #print nbe, cmpname[i], k, pcor, pval
                
                bpval = calc_corr_permtpval(permut_matrix,pcor)
                
                #print nbe, cmpname[i], k, pcor, pval, bpval 
                if bpval <= pvalth:
                    print nbe, cmpname[i], k, pcor, pval, bpval
                    outa[0,0],outa[0,1] = i, k
                    outa[0,2],outa[0,3],outa[0,4], outa[0,5] = pval, pcor, bpval, int(cmpname[i])
                    #print outa
        
                    outa_cor = np.vstack((outa_cor,outa))
                    #print outa_cor
                
                    npairs+=1
            this_iter+=1
    
    print 'The number of accepted correlation is '+ str(npairs)
    
    return outa_cor[1:,]



def calc_correlation_matrix(mat1,mat2, corrtype):
    """
    mat1, mat, are numpy masked arrays, with the same number of columns
    mat1=compounds
    mat2=pathways
    but the different (equal) number of rows, mat1 > mat2.
    Columns correspond to the same variable, ex: sample
    kend, pear: kendaltau or pearson correlation
    """
    m1 = mat1.shape[0]
    m2 = mat2.shape[0]
    parametric = {'Kend':'Kendall','Pear':'Pearson','Sper':'Spearman','Robus':'Roboust'}[corrtype]
    outa_cor=np.empty([mat1.shape[0],mat2.shape[0]])
            
    totc=m1*m2
    this_iter=0

    for i in range(m1):
        nbe =totc-this_iter
        print nbe
        for k in range(m2):
        
            GI50 = mat1[i,:]
            path = mat2[k,:]
                        
            if i < m1:
                pcor, pval = calc_corelation_arr(GI50, path, parametric)
                #print 'These are the correlation and pval '+ str(pcor) + '  ' + str(pval)
                outa_cor[i,k]=pcor
            else:
                break
            
            this_iter+=1
    
    print 'The dimesion of the correlation matrix is  '+ str(outa_cor.shape[0]) +' '+ str(outa_cor.shape[1])
    
    return outa_cor[1:,]


def calc_autocorrelation_matrix(mat1,mat2, corrtype):
    """
    mat1, mat, are numpy masked arrays, with the same number of columns
    mat1=compounds
    mat2=pathways
    but the different (equal) number of rows, mat1 > mat2.
    Columns correspond to the same variable, ex: sample
    kend, pear: kendaltau or pearson correlation
    """
    m1 = mat1.shape[0]
    m2 = mat2.shape[0]
    parametric = {'Kend':'Kendall','Pear':'Pearson','Sper':'Spearman','Robus':'Roboust'}[corrtype]
    outa_cor=np.zeros([mat1.shape[0],mat2.shape[0]])
            
    totc=m1*(m2/2)
    this_iter=0

    for i in range(m1):
        nbe =totc-this_iter
        m2 = np.array(range(m1))
        m2 = list(m2[m2>=i])
        print i, len(m2), nbe
        for k in m2:
            GI50 = mat1[i,:]
            path = mat2[k,:]
                        
            if i < m1:
                pcor, pval = calc_corelation_arr(GI50, path, parametric)
                #print 'These are the correlation and pval '+ str(pcor) + '  ' + str(pval)
                outa_cor[i,k]=pcor
            else:
                break
            
            this_iter+=1
    
    print 'The dimesion of the correlation matrix is  '+ str(outa_cor.shape[0]) +' '+ str(outa_cor.shape[1])
    
    return outa_cor[1:,]



def calc_correlation_mats_mod(mat1,mat2, cmpname, corrtype, pvalth, corth, permutationsFile, pairwise_masking=False):
    """
    mat1, mat, are numpy masked arrays, with the same number of columns
    mat1=compounds
    mat2=pathways
    but the different (equal) number of rows, mat1 > mat2.
    Columns correspond to the same variable, ex: sample
    kend, pear: kendaltau or pearson correlation
    """
    m1 = mat1.shape[0]
    m2 = mat2.shape[0]
    outa_cor=np.empty([1,6])
    parametric = {'Kend':'Kendall','Pear':'Pearson', 'Robs':'Robust'}[corrtype]
    # load libraries for R, just in case the correlation is Robust.
    rstat.load_common_Rlibraries()
    
    
    nper=int(1e5)
    nlabels = int(mat2.shape[1])
    # permuted labels for the GI50 matrix
    #nlabels = int(mat1.shape[1])
    #permtlist = array_permutation(nper,nlabels,arr1=[])
    if not os.path.isfile(permutationsFile):
        print 'Calculating the permutation matrix'
        permut_matrix= corr_permutation_distr(nlabels, mat1, mat2, nper, parametric, pairwise_masking)
        pickledump_matrix_data(permut_matrix, permutationsFile)
    else:
        print 'Loading the permutation matrix ' + permutationsFile
        permut_matrix=pickleload_matrix_data(permutationsFile)
        
    totc=m1*m2
    this_iter=0
    npairs=0
    for i in range(m1):
        nbe =totc-this_iter
        print nbe
        for k in range(m2):
            outa = np.empty([1,6])
            GI50 = mat1[i,:]
            path = mat2[k,:]
            
            if pairwise_masking:
                GI50,path = calc_parwise_masking(GI50,path)

            if i < m1:
                pcor, pval = calc_corelation_arr(GI50, path, parametric)
                #print 'These are the correlation and pval '+ str(pcor) + '  ' + str(pval)       
            else:
                break
            
            if pval <= pvalth and corth <= abs(pcor):
                #print nbe, cmpname[i], k, pcor, pval
                
                bpval = calc_corr_permtpval(permut_matrix,pcor)
                
                #print nbe, cmpname[i], k, pcor, pval, bpval 
                if bpval <= pvalth:
                    print nbe, cmpname[i], k, pcor, pval, bpval
                    outa[0,0],outa[0,1] = i, k
                    outa[0,2],outa[0,3],outa[0,4], outa[0,5] = pval, pcor, bpval, int(cmpname[i])
                    #print outa
        
                    outa_cor = np.vstack((outa_cor,outa))
                    #print outa_cor
                
                    npairs+=1
            this_iter+=1
    
    print 'The number of accepted correlation is '+ str(npairs)
    
    return outa_cor[1:,]


def calc_parwise_masking(arr1,arr2):
    """
    If a value is mask in one of the arrays, it gets mask in the other array also.
    """
    arr1_mask = np.ma.getmask(arr1)
    arr2_mask = np.ma.getmask(arr2)
    common_mask = np.ma.mask_or(arr1_mask, arr2_mask)
    newarr1 = np.ma.array(arr1,mask=common_mask)
    newarr2 = np.ma.array(arr2,mask=common_mask)
    
    return newarr1, newarr2


def corr_permutation_test(arr1, arr2, permlist, parametric, corr):
    """
    Permutation test for the correlation btw pathway activity and GI50
    """

    nper = len(permlist)
    mtrax = np.zeros([nper,2])

    for i,nparr in enumerate(permlist):
        permut_arr2 = arr2[nparr]
        mtrax[i,0], mtrax[i,1]  =  calc_corelation_arr(arr1, permut_arr2, parametric)
        #print 'Permutation  = '+str(i)

    bpval = len(mtrax[abs(mtrax[:,0]) > abs(corr),0])/float(nper)
    
    return bpval


def corr_permutation_distr(nlabels, arra1drug, arra2path, nperm, parametric, pairwise_masking):
    """
    It calculates a permutation distribution for the values of the correlation
    """
    mtrax = np.zeros([nperm,2])
    
    for i in range(nperm):
        
        r1, r2 = np.random.random_integers((arra1drug.shape[0] -1)), np.random.random_integers((arra2path.shape[0]-1))
        # Permute the labels in order to destroy correlation
        permlist = np.random.permutation(nlabels)
        # Get any two arrays from the data matrices 
        GI50 = arra1drug[r1,:]
        path = arra2path[r2,:]
        permut_path = path[permlist]
        #permut_GI50 = GI50[permlist]
        
        if pairwise_masking:
            mGI50, mpath = calc_parwise_masking(GI50, path)
        
        mtrax[i,0], mtrax[i,1]  =  calc_corelation_arr(mGI50, mpath, parametric)
    
    return mtrax

        
def array_permutation(nper,nlabels,arr1=[]):
    """
    Returns generates a permutation list based on an array or its labels
    """
    if arr1:
        permlist1 = [np.random.permutation(arr1) for i in range(nper)]
    else:
        permlist1 = [np.random.permutation(nlabels) for i in range(nper)]
    
    return permlist1


def calc_corr_permtpval(mtrax,corr):
    """
    It calculates p value for the correlation value given for corr.
    """
    
    bpval = len(mtrax[abs(mtrax[:,0]) > abs(corr),0])/float(mtrax.shape[0])
    
    return bpval



def calc_corelation_arr(arr1, arr2, parametric,rlibrary=True):
    """
    Calculates the correlation between two vectors. The vectors are masked arrays
    """
    
    ma1=np.ma.count(arr1, axis=None)
    ma2=np.ma.count(arr2, axis=None)
    # minimum number of observations required to calc correlation
    thr=0.4
    min_obsv = np.ceil(thr*arr2.shape[0])
    # handles what happens when one or both arrays are only masked values
    if ma1 <= min_obsv or ma2 <= min_obsv:
        pcor, pval=0.0,0.0
        print 'Arrays do not have enough values to calculate correlation'
    
    else:
        
        if parametric == 'Pearson':
            pcor, pval = sp.stats.mstats.pearsonr(arr1, arr2) 
            pval = pval.data       
        elif parametric == 'Kendall':
            # If use_missing=True assigns the average value to the missing values, if
            # False it gives them rank 0. Returns
            pcor, pval  = sp.stats.mstats.kendalltau(arr1, arr2, use_ties=True, use_missing=False)
        
        elif parametric == 'Spearman':
            # It computes a Spearman non parametric correlation test
            pcor, pval  = sp.stats.mstats.spearmanr(arr1, arr2, use_ties=True)
            
        elif parametric == 'Robust':
            nparr1 = np.ma.filled(arr1, fill_value=np.nan)
            nparr2 = np.ma.filled(arr2, fill_value=np.nan)
            print np.ma.count(arr1,axis=None), np.ma.count(arr2,axis=None)
            
            arrmat = np.transpose(np.vstack((nparr1,nparr2)))
            
            # Calc correlation using a robust algorithm. rlibrary=True avoids loading the library
            pcor = rstat.calc_robust_correlation(arrmat,rlibrary)
            pval = 0.0
            #Note: there is an error produced by covRob, that do not understand but It is produced when small number of observations
            # to calculate the correlation.
            # try:
            #except rinterface.RRuntimeError:
               
          
    return pcor, pval



def order_arrayMat_bysamples(expMat, samples1, samples2):
    """
    This functions order the columns in expMat (cmpMat) according to the list of samples provided in samples1. 
    It is used to order the compound matrix object in the same way that the path activity matrix is ordered.
    This is necessary in order to compute correlations between pathway activity in sample vs GI50
    """
    
    odsp=[]
    matdim=expMat.expVal.shape

    odexpMat=cpmMatrix(matdim[0], matdim[1])
    spnotfound=[]
    for i, sp in enumerate(samples1):
        if sp in samples2:
            #spind = expMat.sampleInd[sp]
            spind = samples2.index(sp)
            odsp.insert(i,(i, sp))
            odexpMat.expVal[:,i] = expMat.expVal[:,spind]
        else:
            spnotfound.append(sp)
            odexpMat.expVal[:,i] = 0
        
    # 
    print 'This samples were not found '
    print spnotfound
    # Making a copy of the old expMat fields    
    odexpMat.probeID = expMat.probeID
    odexpMat.geneSym = expMat.geneSym
    odexpMat.sample = odsp 
    odexpMat.geneSynd = expMat.geneSynd
    odexpMat.genInd = expMat.genInd 
    
    return odexpMat


def get_arrayMat_spfiltered(expMat, samples1):
    """
    This functions returns an arrayMatrix object (odexpMat) that was filtered according to the list of samples provided in samples1. 
    It is used to for example, extract the NCI60 samples.
    samples1 is supposed a subset of samples2=expMat.samplesInd.keys(), if not that column is filled with zeros.
    """
    
    odsp=[]
    matdim=expMat.expVal.shape

    odexpMat=arrayMatrix(matdim[0], matdim[1])
    spnotfound=[]
    samples2=expMat.sampleInd.keys()
    for i, sp in enumerate(samples1):
        if sp in samples2:
            spind = expMat.sampleInd[sp]
            odsp.insert(i,(i, sp))
            odexpMat.expVal[:,i] = expMat.expVal[:,spind]
        else:
            spnotfound.append(sp)
            odexpMat.expVal[:,i] = 0
        
    # 
    print 'This samples were not found '
    print spnotfound
    # Making a copy of the old expMat fields    
    odexpMat.probeID = expMat.probeID
    odexpMat.geneSym = expMat.geneSym
    odexpMat.sample = odsp 
    odexpMat.geneSynd = expMat.geneSynd
    odexpMat.genInd = expMat.genInd
    odexpMat.create_sample_index()
    
    return odexpMat


            
def compare_gene_lists(list1,list2):
    """
    It compares two list of gene symbols.
    The comparison is made from list1 to list 2
    returns the common and not shared gene lists.
    """
    
    gene1= set(list1)
    gene2=set(list2)
    
    common = gene1.intersection(gene2)
    notshared = gene1.symmetric_difference(gene2)
    
    print 'Common Genes '+str(len(common))
    print common
    
    print 'not shared genes '+str(len(notshared))
    print notshared
    
    return list(common), list(notshared) 
    
    
def write_list_genes(outfile, listoflists):
    """
    Writes to a particular files all lists in the master list.
    """
    
    for thislist in listoflists:
        outfile.write(','.join(map(str,thislist)).replace(',','\t')+'\n')


def write_matrix_vals(outfile, matrix):
    """
    """
    
    for i in range(matrix.shape[0]):
        el = map( str,list(matrix[i,:]) ) 
        outfile.write(",".join(list( el )).replace(",", "\t") + '\n')

     
def normalize_matrix_rows(matdata):
    """
    Normalize the row values of matdata between -1 to 1
    """
    for i in range(matdata.shape[0]):
        mat2 = matdata[i,:]
        matdata[i,:] = 2*(mat2/float((np.max(mat2) - np.min(mat2))))
        
    return matdata


def center_normalize_mat(matdata, row=False):
    """
    center and normalize each row of a particular matrix
    """
    if row:
        #normalize by rows
        ind=1
        exp_mean = np.ma.mean(matdata, axis=ind)
        exp_std = np.ma.std(matdata, axis=ind)
        
        # reshape mean and std vectors to matrices
        exp_mean = np.reshape(exp_mean,(matdata.shape[0],-1))
        exp_mean = np.tile(exp_mean,matdata.shape[1])
        #
        exp_std = np.reshape(exp_std,(matdata.shape[0],-1))
        exp_std = np.tile(exp_std,matdata.shape[1])
        #
        exp_norm = (matdata - exp_mean)/exp_std
        
    else:
        #normalize by columns
        ind=0
        exp_mean = np.ma.mean(matdata, axis=ind)
        # Get the standard deviation of each array-sample. So, mean by column
        exp_std = np.ma.std(matdata, axis=ind)
        # Normalized data
        exp_norm = (matdata - exp_mean)/exp_std
        
        
    return exp_norm




'''
inputfile='/home/oabalbin/projects/networks/output/2010_03_20_21_13/bfrm/2010_03_20_21_13_mF_othersp.txt'
outfile = '/home/oabalbin/projects/networks/output/2010_03_20_21_13/bfrm/2010_03_20_21_13_mF_othersp_norm.txt'

matdata = np.loadtxt(inputfile,dtype=float, delimiter='\t')
for i in range(matdata.shape[0]):
    mat2 = matdata[i,:]
    matdata[i,:] = 2*(mat2/float((np.max(mat2) - np.min(mat2))))
    
write_matrix_vals(open(outfile,'w'), matdata) 
'''
'''
gp = gpm.geneparser()
#gen2 = gp.list_of_names(open('/home/oabalbin/projects/networks/output/svd/classifiers/topgen_cmp9_100.txt'))
#gen1 = gp.list_of_names(open('/home/oabalbin/projects/networks/output/svd/classifiers/topgen_400cmp14_SiRNA.txt')) 
#p2p = gp.list_of_names_in_line(open('/home/oabalbin/projects/networks/input_param/p2p_interactome.db')) 
#golub = gp.list_of_names(open('/home/oabalbin/projects/networks/input_param/Golub_ERG_signature'))
#genes_from_kegg = gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/ERG_regulation_KEGG_hugogenes'))
#genes_from_modules=gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/2010_04_07_09_29_erg_targets_geneInModules.txt'))

gp = gpm.geneparser()
genes_from_kegg = gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/All_ERG_chip'))
genes_from_modules=gp.list_of_names(open('/home/oabalbin/projects/erg_regulation/all_genes_in_ETS_classifier.txt'))


print 'gene lists g1 vs g2'
com, nost = compare_gene_lists(genes_from_kegg,genes_from_modules)
print com
print nost
outfile=open('/home/oabalbin/projects/erg_regulation/genes_modules_chipseq_erg.txt','w')

write_list_genes(outfile,com)


mat = pickleload_matrix_data(open('/home/oabalbin/projects/networks/output/bfrm/2010_02_19_23_20/pathcorr_dump_2010_02_19_23_20.txt'))
mi = mat[mat[:,1]==11,]
write_matrix_vals(open('/home/oabalbin/projects/networks/output/bfrm/2010_02_19_23_20/corrcompound_factor11_.txt','w'), mi)
'''